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Abstract. An improved high order finite difference meth- 
od for low Mach number computational aeroacoustics 
(CAA) is described. The improvements involve the con- 
ditioning of the Euler equations in perturbation form 
to minimize numerical cancellation error, and the use of 
a stable non-dissipative sixth-order central spatial dif- 
ferencing for the interior points and third-order at the 
boundary points. The spatial difference operator satis- 
fies the summation-by-parts property to guarantee strict 
stability for linear hyperbolic systems. Spurious high fre- 
quency oscillations are damped by a third-order charac- 
teristic-based filter. The objective of this paper is to ap- 
ply these improvements in the simulation of sound gen- 
erated by the Kirchhoff vortex. 


1 Introduction 

Owing to the high accuracy requirements in the numeri- 
cal simulation of acoustic waves, efficient high order nu- 
merical methods are most sought after in the emerging 
area of computational aeroacoustics (CAA) [28,30]. It 
has been shown that for appropriate high order meth- 
ods, the number of grid points per wavelength can be 
greatly reduced from that of standard second-order spa- 
tial schemes [4]. Low dispersive fourth-order or higher or- 
der schemes have been shown to be the methods of choice 
for linear or weakly nonlinear aeroacoustics in general ge- 
ometries. Complex and CPU intensive schemes such as 
the fifth or higher-order WENO schemes are generally 
considered as the method of choice if complex nonlinear 
aeroacoustic problems are involved. The present study 
is the first of a series of papers [17,18] in an attempt to 
combine several of the new developments [25,5,27,3,19, 
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23,22,31-33] in efficient, highly parallelizable high order 
non-dissipative spatial schemes with characteristic-based 
filters for CAA that exhibit good long wave propagation 
accuracy for linear and nonlinear problems [33]. These 
papers extend the work of [31-33] for CAA. The goal is 
to propose a scheme that minimizes numerical cancel- 
lation errors, and improves nonlinear stability and ac- 
curacy associated with low Mach number CAA. These 
papers utilize the aforementioned new developments in 
an incremental fashion in order to validate the final ap- 
proach. 

The final form of our scheme consists of two levels. 
From the governing equation level, we condition the Eu- 
ler equations in two steps. The first step is to split the 
inviscid flux derivatives into a conservative and a non- 
conservative portion that satisfies a so-called generalized 
energy estimate [3,23]. This involves the symmetrization 
of the Euler equations via a transformation of variables 
that are func tion s of the physical entropy [7]. This split- 
ting of the flux derivatives, hereafter, is referred to as 
the entropy splitting. The split form of the the Euler 
equations was found to require less numerical dissipation 
than its un-split counterpart in association with non- 
dissipative spatial central schemes [32,33]. Owing to the 
large disparity of acoustic and stagnation quantities in 
low Mach number aeroacoustics, the second step is to re- 
formulate the split Euler equations in perturbation form 
with the new unknowns as the small changes of the con- 
servative variables with respect to their large stagnation 
values [25]. Nonlinearities and the conservative portion 
of the split flux derivatives are retained. This perturba- 
tion form was shown to minimize numerical cancellation 
errors compared to the original conservation laws [25]. 

From the numerical scheme level, a stable sixth-order 
central interior scheme with a third-order boundary scheme 
that satisfies the discrete analogue of the integration-by- 
parts procedure used in the continuous energy estimate 
(summation-by-parts property (SBP) ) is employed [27]. 
The discrete scalar product is based on a diagonal ma- 
trix. If the split form of the inviscid flux derivatives are 
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not used, only linear stability is retained. Characteristic 
and nonreflecting boundary conditions (BCs), if needed, 
are imposed at each time step. To suppress the spuri- 
ous high frequency oscillations associated with central 
schemes, a modified version of the characteristic-based 
filter method of Yee et al. [31] is used. The metric terms 
in the general coordinate transformation are discretized 
by the same difference operator as the flow variables 
leading to freestream preservation (uniform flow conser- 
vation) [29] for the conservative portion of the split equa- 
tions. The time derivative is approximated by a 4-stage 
low-storage second-order explicit Runge-Kutta method 
with careful treatment of the intermediate BC at the 
different stages of the Runge-Kutta method to minimize 
loss of time accuracy [2,4,10]. 

The numerical experiments presented in this paper 
consider only the perturbation form of the Euler equa- 
tions. Numerical results to gain nonlinear stability (and 
further minimize the use of numerical dissipation) via 
the entropy splitting will be presented in [17,18]. 

The numerical method has been applied to the com- 
putation of vortex sound. The prediction of vortex sound 
has been one of the most important goals in computa- 
tional aeroacoustics (CAA), because the noise in turbu- 
lent flow is generated by vortices. To verify the high or- 
der finite difference method for the 2D Euler equations, 
here, we focus on the numerical simulation of a single 
Kirchhoff vortex. The Kirchhoff vortex is an elliptical 
patch of constant vorticity rotating with constant angu- 
lar frequency in irrotational flow. The acoustic pressure 
generated by the Kirchhoff vortex is governed by the 2D 
Helmholtz equation, which can be solved analytically us- 
ing separation of variables [16]. The sound generated by 
the Kirchhoff vortex constitutes a new challenging test 
case for CAA. It allows to check numerical methods for 
the Euler equations on 2D polar grids and to test bound- 
ary conditions at the surface of a sound generator and 
at the farfield. 

The outline of the paper is as follows. The perturba- 
tion formulation of the Euler equations [25] is reviewed 
in Section 2. The summation-by-parts (SBP) finite dif- 
ference operator is reviewed in Section 3. The analytical 
solution for the sound generated by the Kirchhoff vortex 
is described in Section 4. Numerical results are compared 
with the analytical solution in Section 5. 


2 Perturbation Formulation of the Euler 
Equations 

In low Mach number aeroacoustics, the changes in pres- 
sure, density, etc. are much smaller than their reference 
values. For example, the acoustic pressure p' is usually 
many orders of magnitude lower than the stagnation 
pressure po* Computing small differences of large num- 
bers on the computer leads to cancellation. The pertur- 
bation formulation introduced in [25] is used to minimize 
numerical cancellation error for compressible low Mach 
number flow. The Euler equations are expressed in terms 
of the changes of the flow variables with respect to their 
stagnation values. Since the velocity in stagnant flow* is 


zero and the stagnation conditions are constant, the Eu- 
ler equations in perturbation form can be written as 

^ + V(pu)'=0, (1) 

- + v • (pu)'u' + Vp' = 0 , (2) 

at 

^l + V-dpHYu 1 + (pH) 0 u')=0, (3) 


where 


p' = P~ Po, (pu)' = pu, ( pE )' = pE- (pE ) o , 

u' = -1^, p' = ( 7 - 1 )[(pE)' - hpu)' • u'], 

(pHy = ( P Ey+p r . 


Here, p denotes the density, u the velocity, E the total 
energy per unit mass, H the total enthalpy, and 7 = 1.4 
the ratio of specific heats for air at standard conditions. 
The “/” and subscript “0” denote perturbation and stag- 
nation variables, respectively. 

Although the present formulation is mathematically 
identical with the conventional conservative form, dis- 
cretizing e.g. Vp leads to cancellation errors, whereas 
these errors are avoided when discretizing Vp'. In Carte- 
sian coordinates, the perturbed 2D Euler equations can 
be expressed as 


aU ; dF[ 3F' 

dt dx dy 

where 
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Here, u' = u is the x-direction velocity and v' — v is the 
p-direction velocity. 

For the treatment of general geometries, a coordi- 
nate transformation (a:(^,p),p(^,r 7)) is used. The result- 
ing transformed 2D Euler equations are 


dl V dF\ dF! 2 

dt + + drj 

where 


(5) 


U' = 

Fi - 

F' = J-^F; + J~ 1 T] y F' 2 , 


with the Jacobian determinant of the transformation 
J” 1 =2 |||^ - |||| , and the metric terms 

' Sy dr) 1 

J fjx — J r ly ~ d£ * 
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3 Numerical Method 

3.1 Summation-by- Parts (SBP) Operator 

For linear partial differential equations, well-posedness of 
the Cauchy problem or initial-boundary-value problems 
(IBVPs) can be proved by the energy method [12,5]. 
The essential mathematical tool in the energy method 
for continuous problems is integration-by-parts 

(u.Vx) = u(l) T i/(l) — u(0) T v(0) - ( u xy v ). (6) 

Here u and v are differentiable d-dimensional real func- 
tions on [0, 1] and not to confuse with the u and v veloc- 
ities of the 2D Euler equations. The (u,t/) = u T vdx 
is the L 2 scalar product and Hull 2 = (11, u) denote the 
L 2 norm. 

As an example, we consider the scalar linear advec- 
tion equation 

ut + cu x = 0 , 0 < x < 1 , (7) 

u(x,0) = f(x) , 0 < x < 1 , (8) 

u(0,0=fl(0. 0 <f, (9) 

where the wave speed c > 0 is constant. Application of 

the product rule and (6) to (7), we obtain the equalities 

ajIM-.OII 2 = 2(u,u ( ) = — 2c(u, u x ) 

— -c(u 2 (l,t) - u 2 (0, <)) = -cu 2 (l,t) + c</ 2 (t) . (10) 

Note that u T = u for the scalar problem (7). Integra- 
tion over a time interval [0, t] shows that the energy 
l||u(.,0|| 2 can be estimated in terms of the initial con- 
dition (IC) and BCs. Thus, the problem is well-posed. 

For the moment, let’s discretized the computational 
domain [0, 1] by N + 1 grid points xj — jh , j = 0, 1, AT, 
with /i = -E. Denote Vj = Vj(t) and Wj — Wj(t) as the ap- 
proximations of u(xj,t) and ^u(xj,t), respectively, and 
v — [v 0 ,vi, and w = [wo, itfi, u?/v] T ' Kreiss 

and Scherer [13], Strand [27] and Carpenter et al. [1] 
constructed high order difference operators Q for 

w — Qv (11) 


where (C?l 6) u)j = £( 55^+3 “ h v i+z + f u i+i “ \ v i~' + 
is the standard sixth-order central dif- 
ference approximation of the first derivative. The forms 
of the 5x9 matrix D = ( djk ) and matrix H can be 
found in [27,3]. Here H is the diagonal matrix defining 
the norm of the SBP operator. Since the SBP opera- 
tor (13) is based on a diagonal norm, its application to 
multi-dimensions is straightforward. 

To closely maintain the order of accuracy of the scheme 
in curvilinear coordinates, the metric terms are discretized 
by the same difference operators as the flux derivatives 
in (5). We approximate the ^-derivatives ^ by the SBP 

operator Q $ (13) and the 77- derivatives by the Stan- 

Z^X 

dard sixth-order central difference operator . In 3D, 
the Vinokur and Yee [29] treatment of the correspond- 
ing metric terms for freestream preservation is recom- 
mended. 

In order to not destroy the SBP property, there are 
different ways in imposing the physical BCs in conjunc- 
tion with the SBP operator to obtain strict linear stabil- 
ity [1,20,21]. The penalty method called ‘‘simultaneous 
approximation term” (SAT) of Carpenter et al. [1] or 
the projection method of Olsson [20,21] are two popular 
SBP preservation approaches. Either approach yields a 
discrete energy relation similar to the continuous energy 
relation. Nonlinear stability can be achieved by applying 
the boundary schemes to the in-going characteristic vari- 
ables via the entropy splitting form of the inviscid flux 
derivatives. For simplicity, we have implemented the in- 
going Riemann invariants without the SAT or the pro- 
jection operator. We use instead the so-called injection 
method, i.e. by imposing them explicitly (cf. section 5) 
which might destroy the SBP property. 


3.2 Time Integration 

The application of the spatial discretization of the per- 
turbed Euler equations in transformed coordinates (5) 
results in a semi-discrete system of nonlinear ODEs 


such that the summation-by-parts (SBP) property is sat- 
isfied, i.e. 

(u, Qv)h = u N v N -u 0 v 0 - ( Qu,v)h , (12) 

where u,u E M' v+1 . The discrete scalar product and 
norm are defined by 

(u,v)h = hu T Hv IMIa = ( u,u)h , 

where H is a symmetric positive definite (iV-b 1) x (iV+ 1) 
matrix. 

We employ the SBP operator, which is third-order 
accurate near the boundary and compatible with the 
standard sixth-order central difference operator in the 
interior. It was derived by Strand [27] and is of the form 

f ~R Sjt=o djkVk yj — 0, — , 5, 


{Q*v)j = | 


(Qi 6) v)j = 6, (13) 

£*=o d N -j, k VN-k ,j = N -5, 


= R(U) ’ (14) 

where U is the vector of the difference approximations 
U' t and R is the vector of two-dimensional spatial dif- 

dv' t aF' / 

ference operators operating on and with each 

element Rjk = -Jj,k + Qv £2] ■ 

For efficiency, the ODE system is solved by a multi- 
stage method [9] 

U (1) = U n + ^ R(U n ) , 

U {2) = U n + ^ R(U (l) ), 

U (3) = U n + ^ R(U (2) ), 

U n+1 = U n + At R(U {3) ) . (15) 

This time discretization is of 0(At 2 ) for nonlinear prob- 
lems and 0(At 4 ) for linear problems. It ha s the same sta- 
bility domain as the classical fourth-order Runge-Kutta 



method. The CFL conditions for the numerical solution 
of §* + = 0, c = constant , with periodic IC and BC 

are (e.g. [15]): 

a < 2.828 for Q i 2) , a < 2.061 for Qi 4) , a < 1.783 
for where a = ^ is the CFL number and 

denotes the standard central /th-order finite difference 
method. The difference in the phase errors between the 
two Runge-Kutta methods will be addressed in a forth- 
coming paper. 

If Dirichlet BCs are properly treated in time, the sta- 
bility limit for a SBP operator with Q$ in the interior 
applied to an IBVP will in general be the same as for 
Q^x applied to a periodic problem, cf. [4], pp. 202-203. 
However, since the boundary values of time dependent 
Dirichlet BCs lack the ’errors’ expected during the dif- 
ferent stages of the classical Runge-Kutta method, the 
mismatches ruin the normal cancellation of errors to fi- 
nal 4th-order accuracy. Instead, the mismatches lead to 
0{At) and O(Ax) at the boundary and 0(Ax 2 ) globally 
independent of the high order finite difference operator 
used. The problem and remedies are discussed in [2,4, 
10]. Here, we choose the remedy by, at every intermedi- 
ate stages of the Runge-Kutta method, the SBP differ- 
ence operator is also employed at the boundary. After the 
completion of the full step of the Runge-Kutta method, 
the Dirichlet BCs are prescribed. Full accuracy in time 
is achieved using this remedy, but the stability condition 
becomes more restrictive. 


3.3 Characteristic- Based Filter 


For long wave propagation of nonlinear systems, even 
without shock waves and/or steep gradients present, spu- 
rious high frequency oscillations are generated by non- 
dissipative central spatial schemes. To suppress these 
unresolved waves due to the sixth-order spatial interior 
scheme, a modified version of Yee et al. [31] high order 
artificial compressibility method (ACM) filter scheme is 
used. The ACM filter is based on Harten’s ACM [6] 
switch but utilized in a different context. In the Yee et al. 
[31,32] schemes, one time step consists of one step with 
a fourth-order or higher central spatial base scheme. Of- 
ten an entropy split form of the inviscid flux derivative 
is used along with a post processing step, where regions 
of oscillation are detected using the ACM as the sensor, 
and filtered by adding the numerical dissipation portion 
of a shock capturing scheme at these parts of the so- 
lution. The idea of the scheme is to have the spatially 
higher non-dissipative scheme activated at all times and 
to add the full strength, efficient and accurate numeri- 
cal dissipation only at the shock layers, steep gradients 
and spurious oscillation parts. Here we employ a simi- 
lar procedure except we simplified the filter by removing 
the limiter. Since the present application is for low Mach 
number CAA, the limiter which is designed for captur- 
ing discontinuities is not necessary. Another simplifica- 
tion for low Mach number CAA consists in using the 
arithmetic average instead of the Roe average. 



Fig. 1. Kirchhoff vortex. 


At the completion of a full step of the Runge-Kutta 
method, the numerical solution U'”*/ 1 is filtered by 

U" +1 = U'" fc +1 - AtJ jtk [D<V + tf.U'K 1 . (16) 

where we use the third-order operator 

D i V' j k = U',* 

with the difference operator t = aj+i/ 2 ,fc“ a j-i/ 2 ,fc 

For our numerical experiments, the filter coefficient k in 
the range of 0 < k 0.05 exhibits the desired prop- 
erty. The Jacobian matrix of the flux in the ^-direction 

can be diagonalized as . 

The columns of R^ are the right eigenvectors of and 

may be found in [29]. The eigenvalues of define the 

diagonal matrix — diag{u^ — c ^ , + c^) , 

where = uJ _1 £ x 4* vJ~ l £ y , and 
c € =c v /(J- 1 &) 2 + (J- 1 f y ) 2 . 

D n V f jk is defined analogously. Whether a characteristic 
filter instead of a scalar filter is absolutely necessary will 
be addressed in a future paper. 


4 Analytical Solution for Kirchhoff Vortex 
Sound 

The Kirchhoff vortex is an elliptical patch (Fig. 1) with 
semi-major axis a and semi-minor axis b of constant vor- 
ticity V x u = (0,0,cj) t rotating with constant angu- 
lar frequency ft — in irrotational flow [11]. The 

2D flow field constitutes an exact solution of the 2D in- 
compressible Euler equations [14]. The acoustic pressure 
generated by the Kirchhoff vortex is governed by the 2D 
Helmholtz equation, which can be solved analytically us- 
ing separation of variables [16]. 

The normal velocity for an almost circular Kirchhoff 
vortex, i.e. a = R(1 + e), b = R( 1 -c), 0 < e C 1, can 
be approximated by [16] 

u n « 2Reftsin(2(8 - fit)). (17) 

Assuming a harmonic time dependence at the angular 
frequency 2ft for the acoustic pressure 

p'(r,6,t) = p(r, 8)e~ ,2{lt , 
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reduces the wave equation to the Helmholtz equation 
k 2 p -f Ap = 0. 

with wave number k = 2i7/c 0 . Separation of variables 
yields the solution for the Helmholtz equations 

p'{r,e,t) = KiAH^ihryW-W), ( 18 ) 

where 3R(z) denotes the real part of a complex number z. 
H ^ is the Hankel function of 2nd order. The constant 

A = p 0 4Ref2 2 
kH^'ikR) 

is determined by the radial momentum equation 

du r dp f 

Po ~dt ~ ~!fr 

using the normal velocity (17) of the Kirchhoff vortex. 
In [16], the farfield approximation of (18) for kr 1 
is shown to coincide with the farfield approximation de- 
rived with Green’s function [8], pp. 126-128. 

An almost circular Kirchhoff vortex generates a sim- 
ilar sound field as an almost circular rotating imperme- 
able ellipse [24]. 

The exact solution of the 2D linearized Euler equa- 
tions for sound generated by the Kirchhoff vortex is de- 
termined by means of the exact solution of the Helmholtz 
equation (18). The velocity is computed by 

u = V<p, 

where the velocity potential ip is obtained from the rela- 
tion 


Integrating over time and using (18), we get 
1 

tp{r,e,t) = / p'(r,9,T)dr 

Po J 

= -~K(A 

Pq — L J 1 1 

For isentropic flow, the density perturbation is obtained 
from 



5 Numerical examples 

In this section we illustrate the accuracy of the high order 
method by several numerical examples. 


5.1 Rotating Kirchhoff Vortex 

We consider a Kirchhoff vortex with radius R — 2m, 
e = 0.00125, J? = 82.5 j. The stagnation conditions are 

p 0 = 1.3^, c 0 — 330™- Thus, the Helmholtz number 
becomes R = kR — 2 QR/cq — 1. 

A polar grid of mildly stretched near the Kirchhoff 
vortex in the radial direction, and uniform in the circum- 
ferential direction of 129 x 24 is used. The periodic BCs in 
the circumferential direction are implemented by 6 over- 
lapping grid points. The exact solution of the linearized 
Euler equations (cf. Section 4) is prescribed as the IC. 
At the circle r — R and at the farfield r = 64.375m, 
characteristic BCs are imposed after the completion of 
a full step of the Runge-Kutta method. With the den- 
sity, velocity and pressure non-dimensionalized with ref- 
erence quantities po , cq and pqCq, respectively, the Rie- 
mann invariants can be expressed as p' — u n , p f — p\ u t 
and p' + w n . Here, u n is the normal velocity and u t is 
the tangential velocity. The in-going Riemann invariants 
are prescribed using the exact solution of the linearized 
Euler equations, while the outgoing Riemann invariants 
are taken from the numerical solution computed at the 
boundary. For example, at r = R } if c > u n > 0 then 

P Un, = (p U n ) computed > P P — (p P ) exact 0, 
Ut = ( U t )exact and p' + U n - (p ; 4- u n ) exact • If < «n < 

0, the Riemann invariants for the acoustic waves are un- 
changed, while those for the entropy and vorticity waves 

become p* - p' — (p' - P^computed and U t = ( U t ) C omputed ■ 

At the time T ~ 200 with the non-dimensional time 
step At = 0.15, the Kirchhoff vortex has rotated 7.5ra— 
dians. Figures 2 and 3 show the computed solution along 
the positive x-axis. The solution without the filter com- 
pares well with the exact solution of the linearized Euler 
equations, except near x = 5, where high frequency oscil- 
lations are visible. These spurious oscillations are elim- 
inated by applying the characteristic-based filter with 
k — 0.025. The filtered and the exact solutions agree 
up to plotting accuracy. The error in the acoustic pres- 
sure in Fig. 3 illustrates that the high frequency oscilla- 
tions without filter are indeed eliminated by the filter, 
except in the vicinity of the boundaries. With larger k, 
e.g. k = 0.05, the deviations near the boundaries slightly 
increase. It is noted that while the characteristic length 
scale was chosen as L = lm in the computation, we use 
R = 2m in the plots, i.e. x = in the computation, 

w'here x* is the x-coordinate in m, and x = in the 
plots, p in the plots corresponds to the non-dimensional 
acoustic pressure p'. 

5.2 Kirchhoff Vortex Started Instantaneously 

In section 5.1, a Kirchhoff vortex rotating forever was 
considered. In this section, the same problem is com- 
puted with the IC of the Kirchhoff vortex starting in- 
stantaneously from stagnation conditions. The ICs are 
p 1 = u 1 — v l — p ! = 0, except for the circle r — R , where 
the exact solution of the 2D linearized Euler equations 
for sound generated by the Kirchhoff vortex of Section 
4 is prescribed at t = 0. At the circle r ~ /?, we use the 



Acouattc prvMura. H-kR-1, e-0. 00125, rotated angte - 7.5 rad 


Acoustte prawura, H ■ k R « 1, r-O 00125, related angia ■ 7.5 rad 




Fig. 2. Comparison of computed and exact acoustic pressure 
for H = 1, t = 0.00125, rotated angle = 7.5rad. 



Fig. 3. Influence of filter on error of acoustic pressure for 
7€ = 1, c = 0.00125, rotated angle = 7.5rad. 


same characteristic BCs in Section 5.1. At the farfield, 
we now use nonreflecting BCs. The implementation pro- 
cedure is, after the completion of a full step of the Runge- 
Kutta method, the in-going Riemann invariants for stag- 
nation conditions are prescribed. This means that, at the 
subsonic farfield boundary p ' — u n = 0 is imposed. For 
u n < 0, p # - p' = 0 and Ut = 0 are imposed as well. 

The analytical solution of Section 4 is not valid for 
the instantaneously started Kirchhoff vortex. It is only 
valid for a Kirchhoff vortex which has been rotating 
forever. However, since no wave is travelling from the 
farfield towards the Kirchhoff vortex, as long as we have 
stagnation conditions in the farfield, we can assume that 
the analytical solution for the Kirchhoff vortex rotat- 
ing for infinitely long time is valid up to the wavefront 
of the instantaneously started Kirchhoff vortex. If the 
wavefront has left the domain without reflection at the 
farfield, the analytical solutions for the instantaneously 


Fig. 4. Acoustic pressure contours without filter for instan- 
taneously started Kirchhoff vortex with T-L = I, € = 0.00125, 
rotated angle = 7.5rad. 


started Kirchhoff vortex and for the infinitely long ro- 
tating Kirchhoff vortex should agree. 

Figures 4-6 show the effect of the characteristic-based 
filter when the wave front has reached r « 16 (at time 
T = 200, At = 0.15). The numerical solution without fil- 
ter exhibits spurious oscillations whereas the characteristic- 
based filter with k — 0.025 agrees well with the analyti- 
cal solution between r = 1 and r « 14. Their difference 
in accuracy is more apparent from the acoustic pres- 
sure shown in Fig. 6 At x ~ 16, w f e see that in general 
the wavefront cannot match the infinitely long rotating 
Kirchhoff vortex solution, because the instantaneously 
started Kirchhoff vortex has zero acoustic pressure down- 
stream of the wavefront. The discrepancies between the 
numerical solution for the instantaneously started Kirch- 
hoff vortex and the analytical solution for the infinitely 
long rotating Kirchhoff vortex are therefore, have phys- 
ical reasons. 

Next, we consider the same instantaneously started 
Kirchhoff vortex as before, but computes for a longer 
time. After the time reaches 500 (same At ~ 0.15), the 
Kirchhoff vortex has rotated 18.75 radians. Now r , the 
wavefront is at r « 38.5. If the wavefront were reflected 
at the farfield r = 32.1875, we should see it at r ^ 25.875. 
The computed acoustic pressure along the positive z-axis 
showm in Fig. 7 is in excellent agreement with the an- 
alytical solution, if the characteristic-based filter with 
k = 0.025 is used. Without the filter, spurious oscilla- 
tions can be clearly seen. It is interesting to note that 
even without the filter, the wavefront generated at t — 0 
by the Kirchhoff vortex has passed through the farfield 
without visible reflection. The quadrupole structure of 
the acoustic pressure is correctly recovered by the high 
order SBP operator with the characteristic-based filter 
( k = 0.025). See Figs. 8 and 9 for the comparison. 





Fig. 5. Acoustic pressure contours with filter for instanta- 
neously started Kirchhoff vortex with 7i = 1, e — 0.00125, 
rotated angle = 7.5rad. 



Fig. 6. Comparison of acoustic pressure for instantaneously 
started Kirchhoff vortex with 7 { = 1, t = 0.00125, rotated 
angle = 7.5rad. 


6 Conclusions 

A spatially high order summation-by-parts (SBP) differ- 
ence operator [27] has been used to solve the 2D Euler 
equations in perturbation form [25] for the sound gen- 
erated by the Kirchhoff vortex [16]. For the instanta- 
neously started Kirchhoff vortex, characteristic and non- 
reflecting BCs are imposed at the completion of each full 
step of the Runge-Kutta method. Spurious oscillations 
are eliminated by a characteristic-based filter similar to 

[31]- 

In order to gain nonlinear stability for the nonlinear 
Euler equations in the hope of further minimizing the 
use of numerical dissipation, future work includes the 
use of the entropy splitting form of the Euler equations 
before the application of the perturbation form. The 



Fig. 7. Comparison of acoustic pressure for instantaneously 
started Kirchhoff vortex with H = 1, e = 0.00125, rotated 
angle = 18.75rad. 
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Fig. 8. Acoustic pressure contours with filter for instanta- 
neously started Kirchhoff vortex with Ti = 1, e = 0.00125, 
rotated angle = 18.75rad. 

wavelet filter sensors [26] and other possible problem- 
independent coefficients of the characteristic-based filter 
will be sought. 
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